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Tagged-particle motion in glassy systems under shear: 
Comparison of mode coupling theory and Brownian Dynamics 
simulations 
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Abstract. We study the dynamics of a tagged particle in a glassy system under shear. The recently devel- 
oped integration through transients approach based on mode coupling theory, is continued to arrive at the 
equations for the tagged particle correlators and the mean squared displacements. The equations are solved 
numerically for a two dimensional system, including a nonlinear stability analysis of the glass solution, the 
so called /^-analysis. We perform Brownian Dynamics simulations in 2-D and compare with theory. After 
switch on, transient glassy correlation functions show strong fingerprints of the stress overshoot scenario, 
including, additionally to previously studied superexponential decay, a shoulder-like slowing down after 
the overshoot. We also find a new type of Taylor dispersion in glassy states which has intriguing simi- 
larity to the known low density case. The theory qualitatively captures most features of the simulations 
with quantitative deviations concerning the shear induced timescales. We attribute these deviations to an 
underestimation of the overshoot scenario in the theory. 

PACS. 82.70.Dd - 64.70.P- - 05.70.Ln - 83.60.Df 



1 Introduction 

The motion of a tagged particle, expressed e.g. through its 
mean squared displacement (MSD), is a well known and 
very intuitive indicator for the dynamics of a system. For 
a single Brownian particle (dilute limit) under shear, the 
MSD is very anisotropic and shows superdiffusive motion 
for the direction of shear [1], an effect called Taylor dis- 
persion. For the shear pointing in x-direction with shear 
rate 7 and varying in the y-direction, the MSDs in the di- 
lute limit for the different directions read (see the precise 
definitions below), 



l2\(7) 
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Shearing speeds up the random (non-afhne) motion along 
the direction of the flow because fluctuations along the 
gradient (y-) direction let the particle experience varying 
solvent flows. Random displacements along the gradient 



^ Present address; Department of Physics, Massachusetts In- 
stitute of Technology, Cambridge, Massachusetts 02139, USA; 
Electronic address:kruegerm@mit.edu 



direction therefore increase the displacement fluctuations 
in flow direction. At higher densities, the situation is not 
as clear and has been studied extensively in the past few 
years in experiments, simulations and theory (mostly in 
low density expansions |2]). Systems near the glass tran- 
sition have only been studied in experiments and simu- 
lations before [31I11IS1E1IZ1IH] • At high densities, generally, 
the MSDs for the directions perpendicular to the shear 
direction have been found diffusive at long times, with 
diffusivities depending on shear rate in contrast to the sin- 
gle particle case in Eq. ([T]): The shear influence can only 
be transformed to the directions perpendicular to shear 
by particle interactions. In [5], it has been seen that the 
MSD for the a;-direction grows indeed cubically in time, 
for a system near the glass transition. Nevertheless, the 
quantitative relation between the different directions has 
not been demonstrated. 

For a system of non-Brownian particles [5] , where the 
particles attain diffusive motion for the directions perpen- 
dicular to shear only due to interactions, the relations for 
the different directions are similar to Eq. ([T]). In contrast 
to Eq. ((T|) , the shear dependent diffusivities are anisotropic 
in general. 

For super-cooled liquids in general, the dynamics of the 
tagged particle (as visualized by the MSD or the incoher- 
ent density correlation function) has been shown to exhibit 
nontrivial features after switch on of shear, connected to 
the shear stress as function of time plITU] . After switch on. 
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the stress reaches a maximum (sometimes referred to as 
static yield stress), where the glass yields, followed by a 
monotonic decay of the stress down to the stationary value 
giving the 'flow curve'. This scenario, called 'stress over- 
shoot', was shown to be visible in the transient tagged 
particle functions, as the MSD is superdifFusive and the 
density correlation function is superexponential right af- 
ter the stress maximum. 

In this contribution, we study the tagged-particle mo- 
tion close to vitrification including shear-melted glasses. 
Wc focus on the transient dynamics after switching on the 
shear, which we analyze by mode coupling theory and in 
Brownian dynamics simulations. Our paper is composed of 
the following sections. In section [31 we introduce the con- 
sidered system and present the derivation of the equation 
of motion for the incoherent density correlation function 
in section |31 Section H] discusses its numerical solution in 
detail, including a /3-analysis and the discussion of master- 
curves for small shear rates. Section[5]is devoted to derive 
analytic expressions for the MSDs, discussing the Taylor 
dispersion near the glass transition. Numerical results are 
given in section [HI Section [7] closes the theoretical part of 
the paper by discussing the waiting time dependence of 
the MSDs after switch on. 

Finally, we show the results of our simulations in sec- 
tion |S1 which in subsections 18.11 18.21 and 18.31 presents the 
density correlation functions, the focus on the dynamics 
near the critical plateau and the master-curves, respec- 
tively. In these subsections, the glassy transient correlators 
will be shown to have the interesting features of shoulders, 
which we attribute to the slowing down of the system af- 
ter the stress-overshoot. Subsection 18.41 shows the MSDs 
for the different directions, demonstrating the validity of 
the relations connecting the different directions as found 
in section [5] 



2 Microscopic starting point 

We consider a system of N spherical Brownian (bath-) 
particles of diameter d, and the spherical tagged particle 
of diameter dg dispersed in a solvent. The system has vol- 
ume V. The bath particles have bare diffusion constants 
Do, the tagged particle D^. The interparticle force acting 
on particle i (i = 1, . . . ,iV, s) at position r^ is given by 
F,; = —dU{{rj})/dri, where U is the total potential en- 
ergy. We neglect hydrodynamic interactions to keep the 
description as simple as possible. These are also absent in 
our computer simulations to which we will compare the 
results. 

The external driving, viz. the shear, acts on the par- 
ticles via the solvent flow velocity v(r) = 7yx, i.e., the 
flow points in a:-direction and varies in y-direction. 7 is 
the shear rate. The particle distribution function >F{r = 
{ri},t) obeys the Smoluchowski equation [TTlfT^ . 



dt^{r,t)^Q^{r,t), 



N.! 



Q^Q, + SQ^J2^^-i^^-'^^-'^-'^^l (2) 



with K ~ 7xy for the case of simple shear. f2 is called 
the Smoluchowski operator (SO) and it is built up by the 
equilibrium SO, /?e = J2i ^i ' [^i ~ ^i] of ^^e system with- 
out shear and the shear term 6fl = — ^^ di ■ k ■ Vi. We 
introduced dimensionless units, where lengths, energy and 
time are measured in units of d, ksT and d^ /Dq, respec- 
tively. The effect of shear relative to Brownian motion is 
measured by the (bare) Peclet number Pcq = jd^/Do, 
which in these units agrees with the shear rate. 

The formal H-theorem [13] states that the system reaches 
the equilibrium distribution if'e at long times, viz. ile^e = 
0, without shear. Under shear, the system reaches the sta- 
tionary distribution tZ/g with i?^'^ = 0. Ensemble averages 
in equilibrium and in the stationary state are denoted 



(3a) 
(3b) 



dr^ein . . . , 

{...)^^^^Jdrirg{r)..., 



respectively. 



3 Equation of IVIotion for the Transient 
Incoherent Correlator 

The information about the average dynamics of a tagged 
particle is contained completely in the so called incoherent 
density correlator. Under shear, one can define different 
dynamical correlation functions, as discussed in Refs. [TUJ 
|B]. We will start in this section with the transient one, 
for which the external shear is switched on at t = 0. It 
is the general strategy in the MCT-ITT approach (an ex- 
tension of mode coupling theory (MCT) [T3] to sheared 
systems, where ITT stands for 'integration through tran- 
sients' [12]), to start in deriving the transient quantities. 
In the coherent case this is justified by the generalized 
Green Kubo relations for the stress [T5] and the fact that 
the transient correlator can be obtained with the equilib- 
rium structure factor as only input. Here it is a natural 
continuation to derive the equation for the transient in- 
coherent correlator, since we will be able to use many in- 
sights gained from both the coherent and the equilibrium 
case. Furthermore, this approach will lead to the station- 
ary mean square displacements (see section [7]), one of the 
main goals of this contribution, and the transient incoher- 
ent correlator can serve to derive other observables in ITT, 
in the future. The transient incoherent density correlator 
^q(i) (the intermediate scattering function) is defined as 



K(t) 



-?qr 



^^nH^iqit)-: 






(4) 



with the particle position r^. In contrast to the coher- 
ent case, the normalization of the correlator is unity since 
g-iqrsgiqrs _ ^ holds. On tlic right hand side the ad- 
vected wavevector, a specialty of the ITT-approach [TS] 
appears. It reads 



4=1 



q(0 



•ytq^ey 



(5) 
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It appears in Eq. (U) because of translational invariance of 
the infinite system [16] . All wavevectors other then Eq. ([5|) 
lead to zero in Eq. ^ p!^[T5] . Due to this advection, 
the density correlator is, strictly speaking, no autocor- 
relation function for q^ 7^ 0. It can be rewritten using 



quiescent system [T3]. With this, Eq. ^ can be rewritten 
such that the well behaved operator appears 



^-SnH^iqrs _ giq(t)ra 






-snU 



C(i) 



,^nU^-snh^^q 



(6) 



with 



f2l{t)^ 



,''''' Q'{t)Qle-''''\ 



(14) 



(15) 



We see that <?q(i) is an autocorrelation function with re- 
spect to the time evolution of 

Uit) — g^^^tg-*'^^'"* — n'lt+snh -sn'h rn\ 

It is worth noting that, if fi\ and SQ'^ commuted, we would 

have <?q(t) = <?g (t), the equilibrium correlator. This is, 
of course, not the case. The following derivation of the 
equation of motion for <?q(t) is analogous to the coherent 
case [T3] and we will therefore be very brief. 

The time dependence of the evolution operator U{t) 
can be found by differentiation, 

dtU{t) = e"'\n^ - 5n^)e-'"'' = e^^'* QI e""^'*. (8) 

We see that the equilibrium operator appears. To proceed, 
it is reasonable to define a Hermitian operator as was sug- 
gested in Ref. [15] . 



At 7 = 0, i^l{t) is perpendicular to density fluctuations, 
which is not the case for 7^0. The part which is not 
perpendicular can be split off by writing 



nl{t)^nUt) + nUt)- 



(16) 



The first part is perpendicular to density fiuctuations, 
P'^f2Q{t) ~ 0, while the other one is not. The two parts 
read 






-sn'h 



(17a) 
(17b) 



with the function S{t) given by [ISjffT^;^, — ~ ^^ FiVi) 

S{t) - 7 f dt'e-^^''cT,ye^"'''. (18) 



^lit) 



„snH 



n\t 



-snH 



(9) 



with 5fi^ = J2i^i ' '^^ ' i^i + ^i)- ^^^ is i'h*^ adjoined 
of — (5i7t in the equilibrium average. It follows that i7j(<) 
is Hermitian in the equilibrium average, because J7| is 
Hermitian in this average |12j . 



Because of S{t), the second part of f^l{t) couples to den- 
sity fluctuations. 

As is done in the equilibrium case, a reduced time evo- 
lution operator is employed which satisfies 



dtUrit,t')^Urit,t')[2l{t). 



(19) 



{rf2iit)g) = {{e-'^^'^n^le-'^'^'g) 



(10) 



{{Qle~'^''ne~'''''g) = {9*nl{t)iy 



Its formal solution is given in terms of a time ordered 
exponential, where operators are ordered from left to right 
as time increases [T51. 



(11) 

And with f = g the above equation also shows that the 
time dependent eigenvalues of fil{t) are real and negative. 
Because of 



UAt,t') = eJj''''''^'\ 



(20) 



(ere 



Sf2U _ I -iq-r, SnU ^ /p-iq(t)-r, _ 



:<(*), (12) 



f2l (t) has identical matrix elements as the equilibrium op- 
erator for the case of density fiuctuations, only the den- 
sities are replaced by their time dependent analogs as we 
will see when regarding the initial decay rate in Eq. (|25p . 
The equations of motion are derived in the spirit of the 
Zwanzig-Mori projection operator formalism J17j . where 
we use the time dependent single particle density projector 



p'it)-J2skt))Kit) 



(13) 



with complement Q'^t) = l-P'^t). We abbreviate P-'(O) = 
P", the well known single particle projector used for the 



We still need the connection between reduced and full evo- 
lution operators given by 

Uit) = Urit,0) + [ dt'U{t')P' nl{t')Ur{t,t'). (21) 

Jo 

Taking its time derivative leads to the useful operator re- 
lation, 

dtU{t) - Ur{t,Q)nl{t) + u{t)P'nl{t) 

+ / dt'U{t')P'fil{t')Ur{t,t')fil{t). (22) 

Jo 

The equation of motion for the desired correlator now 
follows by sandwiching the expressions above with sin- 
gle particle density fiuctuations e^"^'^' . As already noted, 
the operator fll{t) is not perpendicular to these density 
fiuctuations and the first term on the right hand side does 
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not vanish as it docs at 7 = 0. The equation of motion 
hence contains an extra term, 



equation of motion can then (with the use of the theory 
of Volterra integral equations [21]) be written as 



(23) 



it') 



The extra term A^{t) reads 

Z\^(t) - (e-''i---[/,(t,0)r?t(t)e'^q-.) , (24) 

it vanishes only at time i = and grows to lowest or- 
der like jt. An analogous term appears in the equation of 
motion for the coherent correlator in Ref. [15]. As argued 
there, its appearance is the only disadvantage of this ap- 
proach compared to an earlier one (Ref. [H]). In contrast 
to Ref. [12], the initial decay rate is positive; it is equal to 
the equilibrium initial decay rate for advected wavevectors 
(recall that il'l has negative semi-definite spectrum). 

The positivity of the initial decay rate makes the numeri- 
cal analysis of the equations below more stable. The mem- 
ory function M^{t) contains on the left hand side the well 
behaved operator H^, 

A./^(t,i') = - {gt;f2lit')Urit,t')f2lit)gt^) . (26) 

If we knew an approximation for Mq(t, t') in terms of the 
correlator itself, the equation would be closed apart from 
Z\q(i). But MCT approximations for Eq. ([^ are not de- 
sirable, as was discussed in Refs. [T^[T5] : Approximating 
Mq(i, t') in Eq. (|23]). one would have to be very careful to 
obtain an equation which describes slow dynamics. This 
is not the case for Eq. ([^^ below. Because of this, we 
perform a second projection step following Ref. [15]. To 
decompose the reduced SO appearing in Ur{t,t'), we use 
the projector 



P'it) = gp 



(e-/?I(i)e^) 



'niit), 



(27) 



with complement Q''{t). While P^{t) is, strictly speaking, 
not a projector because it is not Hermitian, it is still idem- 
potent, P*(t)P''(t) = P*(<). It is applied in the following 
way. 



nlit) 



nl{t){Q%t)+P%t)) 



e-/?l(t)e^) 



«^l(i)- (28) 



One can then relate Mf^{t, t') to another memory function, 
m^{t,t'), which is governed by the irreducible operator 

^i (t) [TOll^ . The lengthy calculation which leads to the 
equations below is presented in detail in Ref. |15j : The 



^q(i), (29) 

-{g:^*Qiit'n{t,t')ni{t)g:;). 



with the new memory function 

'i^' ^ r^{t)r^{t 

It is governed by the irreducible operator 



(30) 



(31) 



Eq. (j29p has an extra term compared to the familiar one 
known from quiescent MCT [221123] : The term on the right 
hand side arose from A%{t) in Eq. 



Kit)^iQtu,{t,^)n\{t)ff^^ 



(32) 



It also vanishes at i = and grows in leading order like 
7i. It does hence not influence the fast decay onto the 
plateau for 7 — !• 0. This exact set of equations for the 
incoherent transient density correlator is now suitable for 
approximations in order to get a closed equation for ^q(i). 

The first simplification concerns the source term A%[t) 
arising from the stress expression Z'(i) in Eq. (flS]) . In 
Ref. [T3] it is suggested to set E(€) = in leading ap- 
proximation. This leads immediately to ^q(0 = since 
nlj{t) = follows in Eq. pTbt . and with it Z\^(t) = in 
Eq. ([24]). We note the identity 



e'"'' = e''''\l + U{t)), 



(33) 



and hence e 



snU 



^snu 



with S{t) = 0. Approximating 



S{t) = leads also to a simplification of the memory 
function m^lt, t') because J7j reduces to ^nit). With this, 
the time evolution Ui{t,t') becomes 



U?it,t') = e 



J^Us nl(,s)Q^s) 



(34) 



It is finally in the space perpendicular to density fluctu- 
ations, P'U^{t,t') = = U^{t,t')P''. For the memory 
function follows 



= (p^V)l2te-^^^**'[/.(i,i')e''''*Q^(i)^l£'q( 

'<(t')^lQ^(^')e-^*'t/^(i,t')e'^'*Q^(t)^l4(,) 



(35) 

The allowed insertion of Q^{t') on the left hand side can 
easily be verified; inserting P^{t') at the same position 
leads to zero. For the following mode coupling approx- 
imations, the pair density projector is used, which is as- 
sumed to describe the slow dynamics in the glassy regime. 



Matthias Kriiger et al.: Tagged-particle motion in glassy systems under shear 



In contrast to the coherent case, the pair projector in the 
incoherent case consists of the product of coherent and in- 
coherent fluctuations [13] . This has a physical reason; the 
fluctuating force Qgi^^gi = Qsl^Q-Fse"^'""") on the tagged 
particle depends on the tagged particle and the collective 
dynamics, i.e., the dynamics of the surrounding bath par- 
ticles. Technically this is achieved by the projector. 






gp(t)gk(t))(gp*(t)gk(t) 



(36) 



with gk = J2i e*"^""' the density of the bath particles and 
Sk = {Qi^Q'k)/N the structure factor. Note that in contrast 
to the coherent pair projector, the two densities can be 
distinguished here and the wavevectors are not ordered. 
No counting factor will appear. The memory function p5[) 
is written. 



m^(i,^') 



1 



7rM^it')^^^Q-'it')Piit')' 



-Sf2W 






(37) 



and in accordance with Ref. |15j , the appearing four point 
correlation function is approximated as the product of cor- 
relators with full dynamics, 






NSk^t')'l>Ut')(i-i')'^W'}it-t')Sp-P'S- 



k.k'- 



(38) 



This factorization of the four point function is the major 
approximation in this approach. A similar approximation 
is used also in quiescent MCT [Tl]. The remaining parts 
of the vertex are now found easily, since they are identical 
as in equilibrium using the advected wavevectors instead 
of the time independent ones. The vertex in equilibrium 
reads (we have already inserted the restriction of p = 
k-q), 



Vqk = 



Qi.^g^kQ' OIq'J ^ 



NSk 



= -^kqc|. 



(39) 



where c* = {g't* Qq) /(nSq) is the direct single particle cor- 
relation function [21], n = N/V is the density. The sum 
over bath particles does not contain the tagged particle, 
and we have rf = {Sg — l)/{nSq) if the tagged particle is 
identical to the bath particles. Summarizing, we find the 
following approximate equation of motion for the incoher- 
ent transient density correlator. 



dt<i>^^{t) + r^{t) S^<p^{t) + J* dt' 



,<(t,i')9*"^q(t') =0 



with r^{t) = q^{t) (compare Eq. ^) and 



(40) 



mUt.t') 



k(f) . q(i) k(i') . q(i') 



N ^ q^{t) q^t') 



^^4it)4it')Skit') Kit')-qit')it ~ i')^k(t')(^ - *')• 



Changing the summation index from k to k' = k(t') (and 
immediately renaming the dummy variable fronr k' to k), 
we get 



S(M') = ^E 



k(t - t') ■ q(t) k • q(t') 



q'it) 



q'it') 



n\' 



k(t-t') 



4SkK^^^t')it-t')<Pu{t-t'). (41) 



We see that this final form only depends explicitly on t' 
via q(i) since we can use, e.g. q(i) = q(i')(^ ^ ^') to write 
m^(<,i')=™^(t-)(i-i') with 



<(i 



k 



1 V- k(i - t') ■ q(t - t') k • q 



q'-it-t') 



n\' 



k{t-t 



,flSk1>l_^{t-t')1>^{t-t') 



(42) 



Through the pair density projector, the dynamics of the 
incoherent correlator is coupled to the coherent correlator. 
Eq. (jM)) can therefore only be solved if the corresponding 
equation for the coherent dynamics has been solved before. 
This coupling is physically intuitive, since a (large enough) 
tagged particle can only move if the surrounding particles 
move. There is a certain percolation threshold for the size 
of the tagged particle, below which it is mobile even if 
the bath is arrested |25j . Yet, we will in the numerical 
solutions consider the case where the tagged particle is 
one of the bath particles, (i.e., the tagged particle is much 
larger then the percolation threshold). Then, at 7 = 0, 
the dynamics of the tagged particle follows the dynamics 
of the bath particles [^[Hll^ . i.e., the tagged particle is 
trapped if, and only if, the bath is arrested. 

The memory function (|4ip depends on t' and t~t' . This 
complicates the following analysis because the convolution 
theorem cannot be applied. It probably originates from 
the fact that we investigate the transient regime which 
is not time translationally invariant. An equation for the 
stationary correlator should contain a memory function 
depending ont — t' only. 



4 Results for the transient incoherent 
correlator 

4.1 Numerical details 

Let us turn to the numerical evaluation of Eq. (|40[) which 
we performed in D = 2 dimensions for a system of equal 
sized hard discs {ds = d). The only thermodynamic con- 
trol parameter is the area fraction 77 = ^. 

The solution for D = 3 is as yet numerically too costly 
in computer time and memory. For D = 2, we used a 
spherical grid with 100 points in radial direction, q = 
0.2,0.6, 1.0, .. .,38.8. The angular space was divided in 
96 portions, giving a grid oi 9 = 0.065, 0.13, . . . , 2it. The 
number 96 is often divisible by 2 allowing us to give the 
correlator for angles 9 = 7r/2,7r/4,7r/8 and so on, which 
are the most interesting to be analyzed. Note that this 
grid is different compared to the one used in |27j . where 
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(only) the coherent density correlators were determined. 
While the resulting solutions are very similar, the current 
grid has the advantage that the correlators for constant 
q can be given for all 9, so anisotropy effects can be well 
studied. This is not possible for the Cartesian grid used 
in [57]. On the other hand, the numerical algorithm for 
the spherical grid involves more interpolation procedures, 
since the vector q — k is not on a grid-point. 

From our discretization follows the critical packing of 
77c = 0.6985658 and the exponent parameter A = 0.7155. 
The latter determines all power-law exponents of the the- 
ory. These values differ slightly from the ones found in 
Ref. [ig (77c = 0.696810890 and A = 0.7167) due to the dif- 
ferent discretization of g-space, which is finer in Ref. [55]. 
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4.2 Correlator ^^ 

As noted above, Eqs. (|^0)) and (|¥T|) (together with the 
coherent analogues [l^) show the well known bifurcation 
scenario connected to the glass transition at 77c, separating 
the control parameter region with intrinsically crgodic cor- 
relators from the one where the correlators only because 
of flow decay to zero at long times. 

While this transition is a cooperative effect, i.e., it hap- 
pens for all wavevectors q at the same density, the shape 
of ^q(i) (for both with and without shear) depends on q. 
For densities below the glass transition, i.e. e = ''^''^ < 0, 
the correlator for the system without shear decays to zero 
with time scale Ta, the so called a-relaxation time jl4| . 
The effect of shear does then depend on the dressed Pcclet 
or Weissenbcrg number Pe= JTa- For small shear rates, 
the effect vanishes, 



lim^^^(t)^<?f)(i), 



Taj~iO 



liquid. 



(43) 



This is demonstrated in Fig. [T] where the correlator for 
a liquid state {e = —10"'^) is shown at different shear 
rates. For large Pe, the final decay is dominated by shear, 
and the correlator is anisotropic in q-space, whereas the 
curve with the smallest shear rate shown {Pcq = 10^^) 
is indistinguishable from the equilibrium curve and the 
correlator is isotropic here. 

Above or at the critical density, the correlator of the 
system without shear stays on the plateau characterized 
by the non-crgodicity parameter /|, 



lim<l>f)(t) = /,^>0, 



glass. 



(44) 



At the transition, /^ jumps discontinuously from zero to 
a finite value, given the size of the tagged particle is not 
too close to the percolation threshold [13]. The system 
under shear, however, is always ergodic, since shear melts 
the glass, and <?q(i) decays to zero for any finite 7. Since 
glassy systems are frozen in without shear, the final decay 
from the plateau to zero is governed solely by shear, for 
arbitrarily small 7 — > 0. The dressed Peclet number is 
always infinite because the intrinsic r^ is formally infinite. 



Fig. 1. Transient incoherent density correlator for e — —10^'^ 
(liquid) and qd = 6.6. Shear rates are 7 = 10^" with n = 
2, . . . , 6. For the three largest shear rates, we show the four 
characteristic directions, for the small rates, only q^ = is 
shown for visibility. The curves for the two smallest rates co- 
incide. 




t Dq / d'= 



Fig. 2. Transient incoherent density correlator for e — 10^'' 
(glass) and qd = 6.6. Shear rates are 7 — 10^" with n — 
2, . . . , 9. For 7 = 10~^, we only show the curve for qx — 0, 
together with a fitted compressed exponential with exponent 
/i = 1.05 (dots), see Eq. ^ . 



Fig-IUshows the correlator for a glassy state (e ~ 10^"^) 
at different shear rates. It is seen that the effect of shear, 
and the anisotropy in q-space, prevails up to arbitrary 
small 7. For e > 0, the functions approach a master func- 
tion for 7—^0 and jt = const., which depends only on jt. 
This can be seen in Fig. [2] and will be discussed in more 
detail in sec. 14.41 For the range of shear rates shown in 
Fig. [21 the anisotropy depends hardly on the shear rate, 
probably because even 7 = 10~^ is already quite well de- 
scribed by the 7 — > master function. 



4.3 /3-Analysis 

Further insight into the dynamics near the critical plateau 
can be gained by the so called /3-analysis. It is a non-linear 
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stability analysis of the frozcn-in structure and consists of 
an expansion of the equation of motion, Eq. (j40p , around 
the critical plateau value f^'^ = fq{^ = 0) [H] defined in 
Eq. pi)) . Ref. [22] presents this analysis for the coherent 



transient density correlator <?q(i) 



(PnC 



nU, 



))/{pl 



'•iV-J — \Pq'^ Pq(t)//\PqPq/: 

with yOq = ^^ e*'!""', which can be written near the critical 
plateau as 

<^q(i) = r, + h, {g{t) + C^f-"^-) (i)) + 0{e). (45) 

hq is called the critical amplitude. The dynamics near the 
critical plateau is given by a g- independent isotropic part, 

Q{t), and an anisotropic part, Qq^''°{t). The equation of 
motion of the former is referred to as the /3-equation [29] , 



c(^)(7t)2 + xg^{t) =. ^ / dt'g{t - t')g{t'). (46) 



0.38 



0.36 



0.34 



(3 




0.32 



5 10 15 20 25 30 35 40 
qd 

Fig. 3. Qq''^"''"° {t)/jt at Qx — Qy (identical to ia(|q|)) as 
function of |q|. 



Where e = Ce = C{r] — rjc)/r]c with C ~ 2.1 describes 
the distance from the transition point |16| . For our grid, 
we find A w 0.7155 and c'^' w 3.4, see e.g. Ref. [H] for 
the definitions of these quantities. Note that Eq. (|46)) is 
nonlinear (quadratic) at the critical point. This is ex- 
plained in detail in Ref. |14| . The short time behavior 
of Q{t) must be matched to the short time dynamics of 
the correlator, Q{t — >■ 0) = (to/t)'^, where the match- 
ing time to is determined by the coherent initial decay 
rate. The critical exponent a obeys (with the r'-function) 
A = r2(l - a)/r(l - 2a). From Eq. (gB]), we see that the 
/3-correlator is of order ^/e and |7t|, and we keep our dis- 
cussion to these orders. See Refs. [^[50] for more details 
on the two parameter scaling relation for jt and e. The 
/3-correlator takes for £ > the solution for long times 



g{t:^tb) 



3(7) 



A 



m 



(47) 



Eq. (|47)) describes the initialization of the final shear in- 
duced decay from the plateau to zero. One has i^ = Ve/lil 
for e > and t^ = to for e = 0. The shear independent 
decay from the plateau for the liquid case can be found in 

Refs. [nisiiisa. 

The anisotropic term in Eq. (j45p has been overlooked 
in Ref. |29j . Since the /?-analysis for the incoherent tran- 
sient correlator depends on the coherent one (isotropic and 
anisotropic), we will only discuss the results here. The de- 
tailed derivation of both coherent and incoherent terms 
will be presented in a forthcoming paper. 

Wc consider the case of e > 0, because for e < 0, the 
dynamics is independent of shear for 7 -> and the equi- 
librium discussion is recovered |22] . Expanding the inco- 
herent correlator near the critical plateau (for < £ ^ 1 
and 7i <i; 1), we find that the /3-correlator contains an 
isotropic part given, as in the coherent case, by g{t), as 
well as an anisotropic part Qq'°'"'^'^° (t), 

'?q(i) = /r + K (^(^) + ^^'^'""■""HO) + 0{e). (48) 



The critical amplitude h^t is equal to the one at 7 = [32] . 
The anisotropic term comes from the lowest order terms in 
7t of the memory function mi^{t,0). Here, Qq^''°{t), the 
anisotropic part of the coherent /3-correlator contributes. 
We find that g^'°-"'^^°) (^'^ jg Hi^car in jt and proportional 
to qxqy, 



-t{s,aniso) 



it) 



«(|q|)^7^ 



om' 



(49) 



The term q^Qy represents the expected "quadrupole" -de- 
pendence. For QxQy > 0, the dynamics is slightly slower 
than on average and for q^qy < it is slightly faster, 
i.e., a(|q|) > holds for all |q|. The function a(|q|) in- 
creases slowly with q, see Fig. |31 The maximum value of 
the anisotropic part (on our grid) is at roughly 0.377^. 
Still, it renders the slope of (P^ positive for the region 
qx ~ qy, since the isotropic contribution g(t) is initially 
proportional to (7*)^ (before Eq. (|T7|) holds). 

The fact that the anisotropic part is in lowest order 
proportional to qxqy is not unexpected. There are other 
examples where such a term emerges, e.g. in the distortion 
of the structure factor under shear; it is in linear order in 
shear rate also proportional to qxqy |331l341l55] for liquid 
states. 

In Fig. [4l we present the agreement between the full 
solutions for <P^ and the /3-correlators near the critical 
plateau. The derived /3-correlator is compared to (^q(i) — 
,fq'^)/h^q for different directions of the wavevector q. The 
positive slope of the correlation function for q^ = qy is 
hardly visible as the anisotropy in the /3-process window 
predicted by theory is rather small. We conclude that 
shear flow frees the particle which would be localized in 
the quiescent glass initially in a rather isotropic process. 

Note that in Fig. SI the shape of the isotropic curves 
(solid lines) is independent of |q|; since Qq'°'"'^^°'[t) ~ 
there, giving rise to the well known factorization property. 
The shapes of the anisotropic curves (dotted lines) on the 
other hand do depend on |q|, i.e., the factorization does 
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but it shows that the correlator indeed obeys the scaling 
described above, also for the ease when the memory func- 
tion does not depend on i — i' only. The reason is that 
the adveeted wavevectors causing the deviation from t — t' 
naturally depend on the strain jt. It can be shown that 
the short time solution of Eq. ([501 at e = is given by 
Eqs. dM]) and (gZl), 



<^q+(i^O) = /,^ 



/iqi, 



(51) 



with /iq 



hUi 



a(|q|)^/c) (secEq. ^). 



Fig. 4. The /3-correlator for the incoherent case. We show 
two glassy states (e = 10""^ and e = 0) with 7 = 10~^. The 
wavevector is g = 6.6 in all curves. For both densities we show 
the directions qx = and qy = (solid lines, lying indistin- 
guishable on top of each other) and the isotropic part of the 
/3-correlator (dashed). Only for e = 10""^, we also show the 
directions qx = qy (upper dotted) as well as qx = —qy (lower 
dotted). Inset: Focus on the anisotropy. Shown is the correlator 
as function of angle 6 at time t — 10^ , referenced to the one at 
9 = 0, for the two densities. The line (through the data) shows 
the result from the /3-analysis for e = which is proportional 
to qxqy = q^ cos{9) sin(^). 



not hold. This statement can also be verified by Fig. |31 
The function Qq"'™''°(t)/jt docs depend on |q|. 



The approach to the master function is exemplified in 
Fig. \5\ where the correlators for a glassy state are plot- 
ted on a rescaled time axis. We characterize the master 
functions by fitting to it compressed exponentials of the 
form 



lim <Pi,(t) 

*0,7t=O(l) ^ 



•^rw^/^cxp 



^{t/r^^^^r-^ 



4.4 Qf-master-curves 



For e > and 7 — >■ with jt ~ const., the correlators ap- 
proach scaling functions ^q+ (t) (with? = \ c^^y{^ — \)it -- 

cjt), which depend only on the timeseale set by 7, i.e., 
they are independent of the short time dynamics set by 
Do [12]. The rescaled time i actually corresponds to the 
accumulated strain since switch-on of shear, and the scal- 
ing law for <?q+(f) expresses that the decorrelation is a 
function of the strain only. These functions obey a sealing 
equation, the so called a scaling equation. Its derivation 
(see App. [^ is complicated by the fact that the mem- 
ory function in Eq. ([M)) is not a function oi t — t', but of 
t and t' separately. Because of this, in the equation be- 
low, derivatives with respect to the adveeted wavevectors 
appear (with mq(t) defined in Eq. ([H])), 

at' \dq{t'/d) i(*/'=)^ 'J ^ ^ ' 

(50) 



(52) 

While the resulting value of the fit parameter fq is very 
close to fq , this equality is not enforced by the fitting pro- 
cedure. Both the resulting relaxation timeseale Tq as well 
as the stretching exponent fiq depend on the wavevector 
and the separation parameter e. In Fig. [SJ we show the 
timeseale for q pointing in y direction as function of |q|, 
for both coherent and incoherent correlators at e = 10"'^. 
The fit has been done with the data for 7 = 10^*^. The 
coherent data are included in order to test and verify the 
good agreement to the data from Ref. ^7\ , which were ob- 
tained on a Cartesian grid. The incoherent values of the 
time scale are as expected much smoother as a function 
of q, while for large q, the two cases approach each other. 
This q dependence of the timeseale of the final decay 
is already visible in the /3-eorrelator; Recalling its solution 
for q = qey in Eq. ([T7)) and rewriting Eqs. ([15|) and ([^ as 
_the first order of an exponential decay from the plateau, 
<I>q{t) sa f^cxp{—ihq/fq), wc cxtract the time scale 



r(7) 



for the coherent, and 






"(7) 



^(7) 

•qey 




(53) 



(54) 



The derivatives with respect to the adveeted wavevectors 
complicates also the numerical solution of this equation. 



for the incoherent case. These curves are also shown in 
Fig. [5] We find that the forms ([55)1 and ([M)) indeed de- 
scribe very well the q dependence of the relaxation time 
scale. While the upper equations yield a prefactor of roughly 

— (^ = 0.252. we achieved the best agreement by set- 
ting it to 0.385. This difference is not unexpected since the 
relaxation time scale depends on e, and we are comparing 
the values for e = (Eqs. ([5^ and ([M)) ) to the one at 
e=10-3 (Fig. [6]). 

The relaxation timeseale of the master-curves depends 
also on the direction of q. This dependence is shown in 
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1/4 1/2 3/4 1 

9/71 [Notation: qx=q cos(e), qy=q sin(e)] 



Fig. 5. Transient incoherent density correlator for e = 10^'' 
and qd — 6.6 e,, (upper curves) and qd — 12.6 e^ (lower 
curves). Shear rates are 7 — 10^" with n = 2, . . . , 9 (for 
qd = 6.6 Gy, same data as in Fig. [J]). Here the time ajcis is 
scaled by shear rate to demonstrate the approach to the mas- 
ter function. Dots show fitted compressed exponentials with 
exponents /i = 1.05 (upper) and 1.13 (lower curve). 



Fig. 7. Relaxation time scale of the incoherent master-curve 
as function of angle for e = 10~^. Full symbols show the 



coherent 

incoherent 

Ret. [27] 




Fig. 6. Relaxation time scale of the master-curve for e = 10^ 
and q — qey. The lines show the time scales as estimated 
from Eqs. (|53|l and (|54|l . The inset shows the small g-data for 
the incoherent case in a logarithmic graph, demonstrating the 
divergence with 1/q^. The line shows the slope of -2. 



Fig. [71 where Tq' is plotted versus the angle 9 (defined 
by Qx = qcosO, Qy = qsinO) for various values of q. We 
see that in most cases, a direction between 9 = n/A and 
9 = tt/2 has the largest relaxation time. While the de- 
pendence on q of the relaxation time scale can be well 
understood by the f3 analysis (compare Fig.|6l), this is not 
quite true for the angular depedence: From the finding 

that Qq'"'"^''"' {t)/^t in Eq. (^S)) is proportional to qxqy, 



we would expect that r, 



(7) 



(a + fo sin 6* cos 6') , where 



a and b describe the relative size of isotropic compared 
to anisotropic contributions. This functional form is also 
shown in Fig. [T] We see that the shape of Tq is quite 
different from this naive expectation, at least for small 
wavevectors, while the curve for the largest wave vector 
shown follows this simple form very well. 



timescales for a fit of Eq. (|52|l to the complete relaxation from 
the plateau to zero, i.e., including the regions where the func- 
tions show the shoulder-like deviations. Open symbols show 
the timescales obtained from fitting Eq. H52|) up to 7t = 1 (ex- 
cluding the shoulders). These are not shown for g = 12.6 since 
the two data sets are indistinguishable. 



For small wavevectors, the correlators develop an angle- 
dependent shoulder at long times, and the shape of the 
curves is very different from a stretched exponential. These 
shoulders are an unexpected feature which is also seen in 
our simulations as shown in Sec. [H For the 7 — 10~^- 
curves used to create Fig. [71 these shoulders start to de- 
velop at roughly t = 10^. Fitting the curves up to i = 10^ 
('short fit') yields the timescales shown as open symbols 
in Fig. [71 One sees that these are closer to the functional 
form (a + b sin 9 cos 9). Furthermore, since the difference 
between 'complete fit' and 'short fit' is a measure for 
the shoulder-like deviation from stretched exponentials, 
we note that the development of shoulders is most pro- 
nounced for small q and for the direction near 9 = Stt/S. 
Following this discussion, we show in Fig. [51 the final de- 
cay for all angles of our numerical grid. Shown arc the 
three wavevectors from Fig. [71 and additionally q = 3. 
For q = 'i and q = 6.6, the shoulders are best visible. 
They are present for a small range of angles (compare 
Fig. [71) . We see that the height of the shoulders increases 
with decreasing wavevector. This can be explained by the 
fact that they appear for all q at roughly the same strain 
(7^ « 1) and the curves with large q relax to zero before 
that time. 



5 Mean Squared Displacements 

Knowing the equation for the incoherent density corre- 
lator under shear, we can now deduce from it the ones 
for the mean squared displacement (MSD) of the tagged 
particle for the different spatial directions and show their 
asymptotic solutions for long times. The transient MSDs 
so obtained describe a particle's motion after switching- 
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Fig. 8. Final decay for all angles of our numerical grid (48 
curves for each wavevector). One can clearly see the shoulders 
for small q. We show different shear rates, as labeled, for visi- 
bility. The four directions of Fig. [2] are coded in the same way. 



on of shear at time < = averaged over equilibrium initial 
conditions. 

Before we start, we have to show the connection of 
the density correlator to the MSD, involving coordinates 
a,b,c,d £ {xs,ys,Zs} of the particle at time i or t = 0. 
This MSD has to be formed with the conditional proba- 
bility W2{rt,r'0), that the system is at state-point F at 
time t after it was at state-point I^' at i = [T31IIS]. The 
MSDs wc will be looking for are of the form 



{[ait)+jtbit)-c{0)-jtdm')=^ 

drdr' [a{r)+^tb{r)-c{r')-^td{r')fW2{rt,r'o) 



(55) 

It is a straight forward calculation to show that this mean 
squared displacement is found by taking the limit of small 
q of the corresponding correlator, 

([a(t)+7t&W-c(0)-7td(0)]2) 

2 _ / ^-iq{c+jtd)^0''t^iq{a+jtb)\ 



lim ■ 

9^0 






q- 



(56) 



From this equation, we will be able to derive the desired 
MSDs. This will be done separately for the different direc- 
tions, since the MSDs will be anisotropic, as was already 
seen in the low density case, Eq. ((TJ. 



5.1 Neutral Direction 

The calculation for the neutral direction is in strong anal- 
ogy to the equilibrium case [551[^. Using Eq. ([55)) . we see 
that we have to expand the correlator for q = qe^ pointing 
in z-direction to get 



Sz'{t)^{[z(t)-z{Q)r) 



2 hm 



i-'Z'kW 



(57) 



(5z^(t) is the transient mean squared displacement of the 
particle in ^-direction. Its equation of motion is achieved 
by expanding (j40p to order q^ and identifying the terms 
via (|57p . The equation is then integrated over time to get, 



1 



1 



t' 



T;Sz^{t) + t-- I dt' I dt"ml{t'-t")dt"5z\t")^0, 

with the memory function in the low q limit (see Eqs. 
and (inSD for the definition of F(k,q,i)) 



ml{t) 



k 



/c,A:,F(k,0,t). 



(59) 



Since m'Kt) has only one time-argument, one can rewrite 
the above equation using the standard trick of partial in- 
tegrations and 5z^{t = 0) = 0, 



Sz' 



it) 



dt'ml{t-t')5z^{t') = 2i. 



(60) 



Eq. (|M| now looks similar to the equilibrium case |22) . 
and its schematic version has been studied before |37l 
[SlITU]. The long time limit of the solution corresponds 
to the small-z part of its Laplace transform 5z'^{z) = 
/p dt e~ " ^ 5 z"^ [t) . The convolution theorem can be ap- 
plied. We find for i — > oo. 



lim 5z^{t) 

t—>oo 



2t 



l + m^{z = 0)' 



(61) 



In contrast to the equilibrium case, 'm^{z = 0) is always 
finite under shear and the MSD is always diffusive at 
long times. In the glass, we have liuij-^om'^^z = 0) ex 
\'y\~'^ (compare Eqs. (J53l54p and the a-scaling equation 
in Ref. [29]) leading to the scaling relation at small shear 
rates, 

lim 6 z^{t)^ 2(3 ^\^\t, (62) 

where the coefficient (S^ = (|7|to^(z = 0))^"'^ is asymptot- 
ically independent of shear rate as 7 — ?> 0. We see that 
the long time diffusivity Dz = /3z|7| is then proportional 
to the shear rate and independent of the short time diffu- 
sivity Dq. Shear flow thus enables the particle to diffuse 
also perpendicular to the flow, which highlights that fiow 
melts the glass. The affine average particle motion decor- 
relates the non-crgodic structural relaxation. It becomes 
ergodic in all directions and for all variables that would 
be non-ergodic in the glass. 

The same linear scaling of the diffusion coefficient with 
7/(i^ is also predicted for sheared non-Brownian particles 
[3 , yet the range of shear rates for these predictions is very 
different. The above analysis holds for Pcq ^ 1, while the 
limit of non-Brownian particles is approached for Pcq ^ 1 
[3]. Presumably, also the physical mechanisms differ. For 
Peo <C 1, shear destroys the localization of particles in 
a quiescent glass and causes structural relaxation. The 
relevant length scale is the localization length that can be 
read off from the quiescent MSD and corresponds to the 
Lindemann length at solidification; often it is connected 
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to the picture of 'cages'. For Pcq ^ 1 shear dominates 
over Brownian motion on all length scales except for in a 
narrow boundary layer close to particle contact. 



5.2 Gradient Direction 

The derivation for the gradient direction is similar to the 
neutral direction. The correlator with q pointing in y- 
direction, q = ge^, is expanded, 

6y\t) EE {[y(t) y{0)r) = 2 lim p^. (63) 

The equation of motion follows analogously and reads 

Sy\t) + f dt'mlit - t')5y\t') = 2t, (64) 

Jo 



with the memory function 



;(t) = 5^fc,(i)fc,F(k,o,i). 



(65) 



Note the slight difference in this memory function com- 
pared to the one in Eq. ((59)) : One of the ky is time depen- 
dent. As expected, the long time limit of 6y'^{t) is given 
by 

(66) 



2t 



lim Sy (t) = — 

This leads to a scaling relation similar to Eq. ([5^ for 
glassy states at low shear rates, 



lim Sy^it) = 2l3y\^\t = 2Dl,^h. 



(67) 



absent. It depends on j/(0) and has to be missing in our 
translationally invariant formulation. 

For finite densities, we can not expect the two defini- 
tions in the first line of Eq. (|55)) to still be identical, their 
difference stays in fact unknown. Our approach naturally 
leads to defining the MSD for the x-direction in terms of 
our transient density correlator, 



We have no reason to expect that the coefficients /3j, = 

{\j\my{z ~ 0))^-^ and Pz are equal, i.e., Dy' = /3y|7| 

will take a different number compared to Dz . Indeed, 
these have been found slightly different in simulations [51 
[S]. Otherwise, the qualitative discussion of the physical 
mechanism behind Eq. (p7)) can be taken over from the 
neutral direction. 



5.3 Flow Direction - Glass Taylor Dispersion 

Concerning the MSD in flow direction, we have to note 
that we are seeking space-translational invariant quanti- 
ties. The expression {[x{t) ~ a::(0)]^) is not translationally 
invariant and hence not appropriate (it depends on y(0), 
see Eq. 11])). Quantities which fulfill this invariance are 
(\x(pi - 'jty{t) - x(0)]2) and {[x{t) + jtyiO) - x{0)]^). One 
can show that the two are identical for small densities 

{[x{t) - jty{t) - x(0)]2) = {[x{t) + jty{0) - ^.(O)]^) 

= 2i?ot+^^o7't'. (68) 

Comparing to Eq. ([1]), we see that the drift term y(0)^7^t^ 
stemming from constant motion with velocity y{0)'yt is 



6x\t) = {[x{t) 



^ty{t) - z(0)]2) 



with 



2hm ■ 

1J-+0 



<?9eJ0 - (e-*«^=e^^^*e*«^=e-T'**«^= 



(69) 



(70) 



This definition agrees with the formal one in Eq. (j56p . The 
equation for Sx'^it) can now be gained by expanding the 
equation for the correlator 'Ple^ (i) in q, 



dtSx^t)+ dt'm°{t,t')dt'6x^{t') 



ra.At) 



(71) 



with 



.°(t,i') = 5: [fc. - jtkyit t')] ^^-^F(k,0,i-O- 



ijt'Y 



Because of (compare Eg. ([25 
2 



EhM^2 + 2ft\ 



(72) 



(73) 



we recover the low density limit of Eq. (j68p using rrL^l (t, t') = 
0, as required for non-interacting particles (infinite dilu- 
tion). Because the memory function in Eq. (j7ip is not a 
function of the difference of its arguments only, the analy- 
sis of the leading long time terms of 5x^ for dense systems 
involves a bit more work compared to the other directions 
above, see App. [Bj We find 



,l™^-'W=3 + 3mO(z = 0) 



(74) 



This result deserves some discussion: It can be regarded 
as the Taylor dispersion for Brownian particles in a shear 
melted glass. The MSD in x-dircction grows cubically in 
time as it does for small densities. The intriguing result 
is that the coefficient for the t^ term is connected to the 
long time diffusivity for the y-direction in the same way 
as in the low density limit. This can be further illustrated 
by writing 



2.2 



lim 5x\t) ^ -5y\t)Yt 



■i y ^ ' 



(75) 



which holds identically in the low density limit, Eq. ([T]), 
and was also found in Ref. [S] for non-Brownian particles. 
We see that this relation comes about because for long 
times, 5x^ is governed by m°(t — t'), see Eq. (P5)) . This is 
physically plausible if we recall the reason for the t^-term: 
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If the particle moves in y-direction, it gets a "boost" in x- 
direction due to the shear flow. It is hence not surprising 
that the t^ term is proportional to Dy', but the result 
that the very same relation holds as in the low density 
limit is nontrivial and unexpected. 

Despite the similarities of the glass Taylor dispersion 
and the low density one, there is an important difference: 
In glasses, the long time term in Eq. ([75]) is independent 
of the bare diffusivity Dq (set to unity here) and obeys 
the yield scaling law. 



hm 5x\t) = -M\W^ 

t—^oo o 



(76) 



again, with the same /3y as in Eq. (|67| . It is also possible 
to derive the next order term in 5x'^{t), see again App.lBl 
It reads 

l + mO(z==0) l + mO(z = 0) ' ^ ' 

Such a term proportional to t^ is not present in the low 
density limit, Eq. (|55)) . It comes about because the mem- 
ory function is not a function of t — i'. Recall that we are 
currently calculating the transient MSD. The stationary 
MSD might not have a term of order t^ for t — J- oo. Note 
that the term in Eq. ([77)) is ex 7^t^ for glassy states. 




Fig. 9. Mean squared displacement for the gradient direc- 
tion for different values of the shear rate, 7 = 10^" with 
n — 0, . . . , 6 and e — —10^"^. The curve for 7 = 10^^ can- 
not be distinguished from the equilibrium curve. The dashed 
curve shows the equilibrium MSD from Ref. [28] . 



The scaling relation in glassy states as 7 ^- follows, 

(81) 



lim 5xy{t)^~Py\^\^t'+0{t), 



with (3y as in Eq. ([S7[) . The sign of 6xy{t) depends on the 
sign of 7, which is expected since inverting the direction 
of shearing corresponds to inverting either x ot y. 



5.4 Cross Correlation 

In the system under shear, there is a correlation between 
x and y which is not present without shear, see Eq. ([1]). 
In our translationally invariant formulation, we define it 
the following way 



5xy{t) = {[x{t) - x{0) - ^ty{t)] [y{t) - 2/(0)] 



(78) 



It can be derived considering the correlator for the diago- 
nal direction q(i = 0) = {q, q, 0) leading to 



6xy{t) 



l-<? 



lim ■ 

(j->0 



g(e^+e„ 



M) Sx^t)+Sy^{t) 



(79) 



See App. IB. 2 1 for the derivation of the long time result of 
Eq. ([7^. The leading order of Sxy{t) is proportional to t^ 
as in the low density case. 



lim Sxy{t) 

t—>oc 



7 



l+?7l0(z 



0) 



-Dl^^fh^t. (80) 



The last step followed with the result for the long time 
diffusion in y-direction in Eq. ([SS[) . We see that 6y'^{t) 
and 5xy{t) are related to each other as in the low density 
limit, except for the minus sign. This sign originates from 
our definition in Eq. ([75[) . Note that defining 5xy{t) ~ 
{[x{t) - x(0) + 7^2/(0)] [y{t) - 2/(0)]) instead would yield a 
plus sign in ([SO)) . The simulations described below also 
give this sign difference depending on definition. 



6 Numerical Results for the Mean Squared 
Displacements 

After having solved the equations for the incoherent corre- 
lator ^q(i), we can solve numerically Eqs. ([M]) . ([7T[) and 
([102p . for the MSD in y and x directions as well as the 
cross correlation. In the 2D numerical calculation we can 
of course not discuss the MSD for the neutral direction. 

In Fig. ini we show the MSD for the gradient direc- 
tion for different shear rates in a fluid state (e < 0). As 
was discussed in Sec. 14. 2[ the MSD approaches the equi- 
librium MSD for 7—^0, the curve for 7 = I0~^ cannot be 
distinguished from it. In Fig. [21 we also show the equilib- 
rium MSD for the same e taken from Ref. [5S] . The slight 
disagreement at long times is due to the different grids 
chosen, as discussed above. 

For the glassy state (e > 0), the long time diffusivity 

Dy' (defined below Eq. ([57[)') is governed by shear for 
arbitrarily small shear rates. In the limit of 7 — > 0, the 
scaling law of Eq. ([57[) is approached with j3y approaching 
a constant. Glass curves for e = 10~'^ are shown in Fig.lTUl 

In Fig. [TTj we finally compare the different directions 
and demonstrate the glass Taylor dispersion Eq. ([75[) . We 
see that the MSD for the x-direction cannot be distin- 
guished from the one for the 2/-direction as long as 7^ ^ 1. 
For long times, with 7^ = 0(1), the two functions sepa- 
rate and the one for the cc-directions approaches the long 
time t^ asymptote from Eq. ([75]) . 
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Fig. 10. Mean squared displacement for the gradient di- 
rection for different values of the shear rate, 7 = 10~" with 
n — 1, ... ,9 and e — 10~^. The dashed curve shows the equi- 
hbrium MSD from Ref. m\. 




Fig. 11. Glass Taylor dispersion: The MSDs for differ- 
ent directions for a glassy state, £ = 10~^ and shear rate, 
7 — 10~^. Thin lines show the long times asymptotes, accord- 



ing to Eqs. (pTll . (j76ll and (|8T 



In Fig. [TTJ also the cross correlation —Sxy{t) is shown. 
It is small compared to 6x^ and Sy"^ for 7^ <C 1 and ap- 
proaches the asymptotic law Eq. (|8l]) for 7^ = 0(1). 



7 Stationary Mean Squared Displacements 

In the previous sections, we derived the equations of mo- 
tion for the transient MSDs. It has been noticed that these 
differ from the stationary ones |6] . Before we discuss these 
differences, let us emphasize the similarities between tran- 
sient and stationary MSDs giving rise to the lowest order 
approximation of setting them equal |15| . For long times, 
when the transient MSD has reached its linear (steady) 
dependence on time, it has to follow the long time dif- 
fusivity of the steady state (the system is then obviously 
in the steady state) . In this regime, transient and station- 
ary MSDs must hence approach each other. Consequently, 
long time diffusivities as extracted from transient or sta- 
tionary MSDs must be identical. The ITT approach of 



deriving the transient quantities thus proves very useful 
here: The results in Eqs. dH]), (|66l), dZHl) and dSO]) hold for 
the stationary MSDs as well. 

In Ref. |10| . an approximate relation between station- 
ary and transient MSDs was derived, which builds on the 
waiting time derivative introduced in Ref. |38[ . For direc- 
tions perpendicular to the direction of shear, we found 
for the stationary MSD 6z'^{t) (or 5y1{t)) in terms of the 
transient one introduced in Eqs. ([57)1 and (|63p . 



5z1 



it) 



5z\t) + a-i5z\i)-5zl{t)). (82) 



dt 



5z1{t) denotes the MSD of the quiescent system without 
shear. The pre-factor a is the normalized integrated shear 
modulus. 



■ ds , 



(83) 



\vxy^xy/ 

with shear stress axy ~ — J2i ^iVi- Fo'" the case of hard 
spheres, the integration in (j83p has to be renormalized 
since the initial value {(JxyfJxy) diverges |12lll0j . During 
this renormalization, a free parameter k enters the equa- 
tion for (7, which is independent of shear rate and den- 
sity. The final expression for the stationary MSD for hard 
spheres is hence given by 



JxyC ^xyf 



ds 



d^k klky{~s)ky S'kS'^^ 
(27r)3 fcfc(-s) 



S?. 



-<^k(-.)W- 



(84) 



The last line followed with the ITT expression for the sta- 
tionary stress, see Ref. [K]. We can now evaluate Eq. ([5^ 
with use of Eq. ([M)) . We used k = 9 x 10"^, as estimated 
from comparing with Ref. [10]. Since Eq. ((82|) holds only 
for directions perpendicular to the shear direction, we can 
only apply it to the y-direction in our 2-dimensional nu- 
merical analysis. This is shown in Fig.[T2]for a glassy state 
at different shear rates. For the five shear rates given, 
we have 7CT = 0.038 ± 0.001 (increasing with shear rate). 
These values compare well to the value of 70- = 0.04 used 
to fit the Brownian dynamics simulation data in Ref. |10j . 
We see, that the difference between transient and station- 
ary correlators is most pronounced at intermediate times, 
whereas for short and long times, the two functions coin- 
cide. 



8 Comparison to simulations 

In this section, we will compare our theoretical MCT-ITT 
results to our simulations. Since MCT-ITT is a quantita- 
tive theory, there is no fit parameter to be adjusted. As 
will be illustrated, the simulation results show many unex- 
pected features, some of which are captured qualitatively, 
but not quantitatively by MCT-ITT. The quantitative dis- 
agreement can be mostly explained by the underestima- 
tion of the overshoot scenario of the stress after switch on 
[TUIIS]. as will be discussed in detail. 
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Fig. 12. Comparison of stationary (dashed) and transient 
(solid) MSDs for the gradient direction for a glassy state, 
e = 10"^, and shear rate, 7 = 10"" with n = 3, . . . , 8. The 
necessity to take derivatives in Eq. (I82|l introduces small wig- 
gles. 



The simulation algorithm is an event driven algorithm 
which describes the dynamics of hard particles, i.e. hard 
discs in the 2D case considered. It has been described in 
detail for 3D in Ref. [3^ and its adaptation to two di- 
mensions can be found in Ref. ^7\. We consider a binary 
mixture of hard discs with the diameters oi di = d and 
db = lAd with equal particle number concentrations and 
a total amount oi N = Ni+ Ni, ^ 1000 hard discs in a 2D 
simulation box of volume V. Thus the packing fraction is 
given by 7] = T^-^{dl + df). For this system, we find the 
glass transition point to be at roughly rj'^ sa 0.7948 |1D] . 

Simulations have been performed at packing fractions 
of ?7 = 0.79 (liquid) and 77 = 0.81 (glass), discretizing the 
time in steps of r^ = 5 x 10~'^ in units of d^ /Dq. 

In the liquid [r] = 0.79), for the transient correlation 
functions, 300 independent initial configurations were pre- 
pared and equilibrated for a time r = 10^ (using, only for 
equilibration, a Newtonian dynamics algorithm) which is 
large compared to the a-relaxation time of Tq sa 10'^, given 
in units of d^ /Dq. For the stationary correlation functions 
600 — 800 independent initial configurations (depending 
on the shear rate) were prepared. According to the find- 
ings in Ref. |10j . stationarity was assumed after 7^.^, > 1, 
where t^ is the waiting time after switch on of shear. 

Above the glass transition density, the preparation of 
transient correlators is nontrivial because the system with- 
out shear equilibrates very slowly. So for the glassy sys- 
tems (r/ = 0.81) we prepared 150 independent equilibrated 
sets for the transient correlation functions. Equilibration 
was achieved by waiting for a period t = 1.6 x 10^, cor- 
responding to an average displacement of the particles of 
half their diameters. After this waiting time, the corre- 
lation functions are independent of waiting time. We es- 
timated the a-rclaxation time to be roughly 20% of our 
equilibration time. This is large compared to the time win- 
dow examined in the following, so that this density can be 
regarded glassy in our simulation time window. 



For the glassy stationary correlation functions, 150 — 
300 independent initial configurations (depending on shear 
rate) were prepared. Again stationarity was assumed after 

it^ > 1. 



8.1 Correlators 

Let ri- [tk] be the position of the i-th particle in the j-th 
of a total M simulation sets for a given time tk ■ Then the 
correlators at time tk for waiting time t^ are calculated 
via 



M 



exp 



i{q-ri^{tk +tw) 



j = l \ 1=1 L 

-qxitkVi, {tk + t^) - q- Ti^ (t^)) 



(85) 



where external shear is switched on at the time origin, so 
that ttu = corresponds to the transient correlator. 

In Fig. [131 we show the transient correlation function 
for the liquid {rj ~ 0.79) at different shear rates. We see 
the close analogy to the theoretical curves in Fig. [TJ For 
large dressed Peclet numbers Pe= JTa, the final decay is 
dominated by shear and the curves arc anisotropic. As 
in Fig. [H the direction q^ ~ is (slightly) slower than 
the qy ~ direction which is slower than the q^ = —qy 
direction for all shear rates. The correlator in q^ = qy 
direction shows a strong shear rate dependent behavior 
in its relaxation time: It decays as fast as the one for the 
q^ = —qy direction for small Fcq numbers, but exhibits the 
slowest relaxation time for large Pcq numbers. The plateau 
in Fig. [13] is lower compared to Fig. [1] which we attribute 
to the bidispersity of the simulations. MCT calculations 
for binary hard discs in two dimensions for the quiescent 
system yield a plateau value of /|''^ ~ 0.61 for qd = 6.5 
while the simulation yields /^'"^ « 0.58. [H]. 

Fig. [HI shows the same functions for the glassy density 
(?/ = 0.81). These curves are in analogy to Fig. [51 Addi- 
tionally to the discussion of the liquid curves, we observe 
the emergence of shoulders for the smallest shear rates: 
For the directions qx = and qy = 0, the correlators dras- 
tically slow down at the end of the final relaxation process. 
We attribute this slowing down to the slowing down of the 
system after the stress overshoot, see below, and Refs. [51 
[To] for the discussion of the shear stress overshoot sce- 
nario. 

Remarkably, we observe these shoulders in MCT-ITT, 
compare Figs. [3 and [1 Figs.[71and[Hlalso show that MCT- 
ITT predicts them most pronounced at a direction be- 
tween qx — qy and qy = 0, this is why they are not as 
clearly seen in Fig. [51 

FigurellHalso presents one MCT-ITT curve for roughly 
the same parameters as the slowest of the simulation curves 
for a quantitative comparison. The small difference in plateau 
heights is as expected (we are comparing simulations for 
a binary mixture to theory for a mono-disperse system). 
Apart from that, the time scale of the initial deviation 
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Fig. 13. Transient incoherent density correlators for rj = 0.79 
(liquid) and \q\d — 6.5 for different directions as labeled. Peo 
numbers are 2x 10°, 6 x 10"\2x 10~\6x 10"^, 2 x 10"^6x 



10"^ 2 
right. 



X 10"^ 6 X 10""*, 2 X 10"^ and Pe, 







from left to 



from the glassy plateau for the Poq = 10~^ curve agrees 
well with that of the Peg = 2 x lO"** simulation curves, 
i.e., MCT-ITT differs at most by a factor of five in shear 
rate. But for larger times, the simulation curves are much 
steeper (compare the compressing exponents in Fig. [T8l 
below) compared to the theory. We attribute this effect 
of large compressing exponents to the stress overshoot 
scenario after switch-on. While the theory curves quali- 
tatively capture this compressing effect (the exponents in 
Fig. [15] are greater than unity), it quantitatively under- 
estimates it. E.g., the memory function in Eq. (|40| does 
become negative for certain parameters (leading to slightly 
negative correlators at long times in Figs. [2]and[T]), but 
the effect is much smaller compared to the simulation. As 
is seen in Fig. [T31 the described underestimation leads to 
a larger deviation of the curves at long times (compare 
also Fig. [T7] below). 

While the stress after switch on for the glassy state 
will be presented elsewhere |42| . in order to underpin the 
conclusions of this section, we marked two characteristic 
times in Fig. \T^ namely the stress (overshoot) maximum 
as well as the time where it has approached its (lower) final 
value. First, we see that the MCT-ITT curve indeed starts 
deviating from the simulation curves at roughly the time, 
where the stress is at its maximum, underestimating the 
successive fast decay. Second, the shoulders indeed start 
emerging when the stress has reached its final value, and 
the dynamics seems to slow down drastically. 



8.2 /3-regime 

Following the discussion in Sect. 14. 3[ it is interesting to 
compare simulation and theory in the /3-regime, where the 
correlators are near the glassy plateau. For simplicity, we 
consider only a single wavevector, focussing on quantita- 
tive comparison rather than testing factorization proper- 
ties. For this test, it is sufficient to regard ^q(i) — /<j, 
where the plateau value fq was chosen appropriately in 



Fig. 14. Transient incoherent density correlators for rj — 0.81 
(glass) and \q\d — 6.5 for different directions as labeled. Peo 
numbers are 2x 10'',6x 10"\2x 10"\6x 10"^, 2 x 10"^, 6 x 
10"^ 2 X 10"^6 X 10"'', 2 X 10"" and Peo = from left to 
right. For Peo = 2 x 10"^ and 2 x 10"'', arrows mark the 
stress-maximum (pointing up) and the time when the stress heis 
reached its final value (pointing down); from Ref. [32] • Dashed 
line is the MCT-ITT result for q = 6.6, q^ ^ 0, £ ^ 10"^ and 
Peo = 10"^ 




t Dp/d 



Fig. 15. The incoherent transient correlator near the plateau 
for 77 — 0.81, Peo = 2 x 10"'* and \q\d = 6.5 (simulations, 
datapoints). Lines show theory curves, with e = 10""^, \q\d — 
6.6, and Peo = 2 x 10"'' (right set) and Peo = 10"^ (left set). 



simulation and theory. Dividing by h^ (compare Fig. [4]) 
is only necessary when testing factorization properties or 
comparing to the /3 correlator, here it would only lead to 
a stretching of the y-axis. Fig. [15] shows the curves for 
q — 6.6, restating that the time scale for the initial decay 
from the plateau is quantitatively described by MCT-ITT 
with a multiplicative error between 1 and 5. Additionally, 
wc observe that the initial anisotropy predicted by MCT- 
ITT ((X qxqyi compare Fig. [4]) is within the statistical noise 
of the simulations and a detailed comparison has to be left 
for future work. 
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Fig. 16. Transient incoherent density correlators for rj = 0.81 
(glass) and \q\d = 6.5 (upper curves) and \q\d — 16.22 (lower 
curves), for the direction q^ = Qy, rescaled by jt. The correla- 
tors approach master functions for Peo — >■ 0, which have been 
fitted by Eq. ((52} (dashed lines). 



8.3 a Master Curves 



As discussed in section 14.41 MCT-ITT predicts the ap- 
proach of a master-curve for small shear rates in glassy 
states (or more precisely in states with Pcq ^ 1 and Pe^ 
1), which depend on time only via accumulated strain 'jt. 
In Fig. 1161 we demonstrate that the simulation curves in- 
deed approach a master function for the system at ?y ~ 
0.81, exemplarily for the direction Qx = Qy] similar behav- 
ior for the other directions can be observed in Fig. HM 

When comparing the properties of the master-curves 
in detail, we first have to note that in both theory and 
simulation, only parts of the a-process can be well fitted 
by a compressed exponential, Eq. ([5^ . On the theoretical 
side, the direction g^; = is an exception, since it can be 
well fitted by a compressed exponential throughout the 
a-process, compare Fig. [51 and we will use this direction 
for the comparison. On the simulation side, the curves are 
almost isotropic up to jt = 0.1 and, e.g., the qx = Qy di- 
rection can be well described by a compressed exponential 
except for the very last part becoming negative and then 
oscillating back. These fits describe all directions up to 
the point were the shoulders emerge. 

In Fig. [T71 we show the comparison of the relaxation 
timescale obtained from this fitting procedure. The the- 
oretical curve is identical to the one in Fig. [51 For the 
simulations, the fit was done with our smallest shear rate 
(Peo = 2 X 10"'*). We see that, while the overall shape of 
the time scale as function of q is the same in theory and 
simulation, the theory overestimates the relaxation time 
by about a factor of 70 for the smallest q. For large q, the 
agreement is much better as there is roughly a factor of 5 
difference. For q = 6.6, the difference is roughly a factor 
20, i.e., more than 4 times larger then the deviation for 
the initial decay from the plateau (compare Figs. [U and 
fTS)) . This additional factor can hence be attributed to the 
underestimation of the compression in the final decay by 
theory. 
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Fig. 

(?7 = 0.81) as function of wavevector. See main text for 
cussion of the differences. 
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Fig. 18. The compressing exponent /iq of the master function 
{rj — 0.81) as function of wavevector. 



This line of argument carries over to Fig. [181 where 
we compare the stretching exponent (in our case rather a 
compressing exponent) for the master functions. As dis- 
cussed before, the fact that the exponent is larger than 
unity can be interpreted as a signature of non-stationarity, 
as it seems to only appear for transient quantities |10l l6]. 
While MCT-ITT correctly captures this nontrivial feature 
on a qualitative basis, the exponent is much larger in the 
simulations. We additionally see that the exponent in the 
simulations has a maximum as function of q, which can 
again be understood as a consequence of the stress over- 
shoot: For large q, the functions have relaxed to zero before 
the overshoot sets in (compare the timescales in Fig. fTT)) . 
hence they do not feel the effect of overshoot and are less 
compressed. As mentioned above, this effect seems not 
to be captured by our theory, as the exponents increase 
steadily with q. Further evidence is given in Fig. ll7l where 
theory and simulation approach each other for large q, 
where the overshoot effect plays no role. We also want to 
emphasize that the direction qx = has for most cases the 
least steep curves, for other directions, we observe expo- 
nents as large as nearly 2 in theory. 
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Fig. 19. Transient incoherent correlators at 77 = 0.81 (glass) 
for the directions as labeled, for wavevectors \q\d — 2.99, 
\q\d = 6.63, \q\d = 8.96 and \q\d = 12.6 from right to 
left. All correlators are at the lowest shear rate accessible 
Peo — 2 X 10"*, but shifted in time by a factor 7 — 10°, 



10' 



10 and 10 from left to right for visibility. 



Fig. 20. Comparison of stationary (dashed) and transient 
(solid) MSDs for the gradient direction at 77 = 0.79. Peo num- 
bers are 2 x 10°, 6 x 10"\ 2 x 10"\ 6 x 10~^ 2 x 10"^ 6 x 10"^ 
2 X 10"^, 6 X 10"'*,2 X 10"" from left to right. Thin dashed 
lines show fits to the asymptote 2Dy t. 



Finally, Fig. \T9\ shows the master-curve for rj = 0.81 
(obtained for Peo = 2 x 10"^) for different wave vectors 
(compare Fig. [5]), allowing to study the wave vector de- 
pendence of the shoulders. For small q, the shoulders take 
up about 40% of the a process, while for larger q, they 
only emerge during the last 20% of the final relaxation 
(estimated from comparing the height of the point, where 
the correlators start to split, to the plateau height). As in 
theory (Fig. [8]), we observe that the shoulders are a co- 
operative effect, happening for all wavevectors roughly at 
the same time -yt ~ 0.07 (the time corresponding to the 
slowing down of the system after stress relaxation), and 
we expect that the effect vanishes for large q (where the 
correlators are zero when this happens). Comparing with 
Fig. m we note that MCT-ITT slightly overestimates the 
shoulders for small q. 



8.4 Mean Squared Displacements 

Let us finally discuss the mean squared displacements. 
Given the definitions 6xi^ (ife, <«,) = x^. {tk+tu,)-'jtkyi^ {tk + 
tw) - Xi^ (i^) and Syt. {tk, t^) = Vi, {tk +tw) - yi- {t^) for 
the displacement at time-step tk for a particle i in x,y- 
direction for the j'-th simulation run, we define the mean 
squared displacements in a similar manner as in Eqs. (|63|) . 
(|M|) and ((75)) . Again 1^=0 and '^t^, = 1 for transient and 
stationary cases, respectively, 

M N 

Sy'^{tk,tw) = jijT^^Syi^{tk,t^)Sy,^{tk,tuj), (86) 



j=i 1=1 

M N 



Sx'^{tk,t^) = ■T-rir;^'^Sxi^{tk,t^)Sx^.{tk,U,), (87) 



j=i 1=1 

M N 



Sxy{tk,t.uj) = — — ^ ^ 6xi. {tk,tn])Syi. {tk, t^ 



3=1 1=1 



The difference between transient and stationary curves 
have been discussed in Ref. [10] and demonstrated in sim- 
ulations for the liquid case. Here we show these curves 
again for completeness in Fig. ^U\\ . We see that station- 
ary and transient correlators approach each other for long 
times. This supports the argument given in Sec. [7] that 
our result (j66p indeed gives the long time diffusivity of 
the steady state. 

Fig. [5T] shows the same plot for the glassy case which 
was not presented in Ref. [T^. Here we can test another 
prediction from the theory; The difference between sta- 
tionary and transient curves prevails up to arbitrarily small 
shear rates (as long as Pc3> 1 holds). This is a nontriv- 
ial statement in agreement with Fig. 1121 There is as yet 
a qualitative difference between theory and simulations 
concerning the transient curves. The simulations show su- 
perdiffusive behavior connected to the stress overshoot 
[TUIIS] which is underestimated in the theory (Figs. [TUl 
and [12]), as discussed before. One possibility for 5y^ in 
Eq. (j64p to show superdiffusive behaviour is a negative 
memory function at long times [5] . But the memory func- 
tion in Eq. (|64p numerically turns out to be positive. There 
is no mathematical reason for this positivity inherent in 
the structure of our equations, and indeed it seems sim- 
ple coincidence that m°(i) is positive for all t: Changing 
its structure slightly can lead to negative values for long 
times and yield superdiffusive motion. 

As was the case for the relaxation time scales in Fig.lTTl 
this underestimation of our theory gives rise to rather 
large deviations of the long time diffusivities from the sim- 
ulation values. In Fig. [^T] we additionally show the the- 
oretical transient curve for Fcq ~ 10"^ demonstrating a 
scenario equivalent to Fig. [T3] While there is in principal 
no free parameter in our theory, we multiply both axis of 
the theoretical data by a factor of 1.22. This factor sets the 
plateau values equal, which are naturally slightly different 
in the binary mixture compared to our mono-disperse the- 
ory. It has no effect on the timescales of the curves which 
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Fig. 21. Stationary (dashed) and transient (solid) MSDs for 
the gradient direction at 77 = 0.81 (glass). Peo numbers are 
2 X 10°, 6 X 10"\2 X 10~S6 x 10"^2 x 10~^6 x 10~^ 2 x 
10"^, 6x10"", 2x10"'' from left to right. The dotted (magenta) 
curve is the theoretical result for Peo — 10^'' and e — 10^^. 
Thin dashed lines show fits to the asymptote 2Dy t. 



Fig. 23. Stationary (dashed) and transient (solid) MSDs for 
the flow direction at 77 = 0.81 (glass). Peo numbers are 2 x 
10°, 6 X 10"\ 2 x 10"\ 6 X 10"^ 2 X 10"^ 6 X 10"^ 2 X 10"^, 6 X 
lO""*, 2 X 10"* from left to right. The asymptotes 2/3D^'''f^7^, 
using the same Dy as in Fig. [2T] for the corresponding shear 
rates, are shown as thin dashed lines. 




tDM" 



Fig. 22. Stationary (dashed) and transient (solid) MSDs for 
the flow direction at 77 = 0.79 (liquid). Peo numbers are 2 x 
10°, 6 X 10"\ 2 X 10"\ 6 X 10"^ , 2 x 10"^ 6 x 10"^ 2 x 10"^ , 6 x 
10"*, 2 X 10"* from left to right. The asymptotes 2/3D^'^'f^7^, 
using the same Dy as in Fig. [20] for the corresponding shear 
rates, are shown as thin dashed lines. 



we want to discuss here: We see that the theory curve for 
Peo = 10"'^ leaves the glassy plateau at the same time as 
the simulation curve for Pcq = 2 x 10"^, in agreement to 
what we observed in Fig [141 Again, up to this point, the- 
ory and simulation agree up to a factor of less than 5 in 
Fcq. Going to larger times, the memory function rriyit) in 
Eq. (|M| misses to become negative and the theory curves 
arc not steep enough. As seen in the plot, the long time dif- 
fusivity differs then by roughly a factor 55 (hence roughly 
10 times more then initially). Regarding again the result 
for the long time diffusivity in Eq. ([66)) . one sees that a 
negative part in my{t) could possibly render my{z = 0) 
much smaller giving larger values for the diffusivities. 

There are three more things to test in our simula- 
tions regarding the flow direction presented in Figs. [35] 
(liquid) and [23] (glass). First, we note that the simula- 



tions indeed show the glass Taylor dispersion, the MSDs 
grow proportional to t^ for long times. Second, transient 
and stationary curves also merge for the flow direction at 
long times, i.e., the expression (|7^ holds for both tran- 
sient and stationary curves as argued in Sec. [7] Third, 
our simulations indeed confirm the nontrivial statement 
of Eq. (|75|) for both hquid and glass: The t^ term is con- 
nected to the diffusivity for the y direction as in the low 
density case, Eq. ([T]). In 3D systems, where diffusivities 
are slightly anisotropic for the two directions perpendicu- 
lar to the flow [B] , we predict that the t^-term is connected 
to the gradient direction rather than the neutral direction, 
as expressed by Eq. (|75|) . 

Inspecting the cross correlator Sxy{t,tw), as shown in 
Figs. [24] and [25] for liquid and glass, respectively, we can 
confirm further predictions of the theory. Again for long 
times transient and stationary functions coincide as ex- 
pected from section [7] Furthermore the connection be- 
tween shear and gradient directions can be seen, as the 
long time asymptote (shown as blue lines) Dy jt^ uses 

the same Dy as in Figs. [20] and [21] for the correspond- 
ing shear rates. This confirms the theoretical prediction 
expressed in Eq. (|M| . 

As a further and more sensitive test of the scaling prop- 
erty in Eqs. ([75]) and §7^, the quantities 3Sx'^{t, t^)/{2'j^t^) 
and Sy'^{t,t^)/{2tj) are shown in Fig. [^S] for the tran- 
sient and stationary curves. For long times the so defined 
mean-squared displacements for shear- and gradient di- 
rection collapse on the constant (3y = Dy j^ for the two 
largest shear rates as expected, while the two lowest shear 
rates already show the right trend, presumably reaching 
their asymptote outside the window accessible in the sim- 
ulation. The inset of Fig. [25] magnifies the gradient di- 
rection for both transient and stationary curves, where 
the superdiffusive regime of the transient mean-squared 
displacement expresses itself by a dip before reaching the 
long time asymptote. For decreasing shear rate, the curves 
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Fig. 24. Stationary (dashed) and transient (solid) cross corre- 
lators at 77 = 0.79 (liquid) . Peo numbers are 2 x lO" ,6x10" 
10~\ 6 X 10"^ 2 X 10"^ 6 X 10"^ 2 X 10" " 



2x 
6x 10"", 2 X lO""* 
^ using the same 



from left to right. The asymptotes Dy jt 

values of Dy as in Fig. [20] for the respective shear rates, are 

shown as thin dashed lines. 




Fig. 25. Stationary (dashed) and transient (solid) cross corre- 
lators at r; = 0.81 (glass). Peo numbers are 2x 10°, 6 x 10^^, 2x 
10"\ 6 X 10"^ 2 X 10"^ 6 X 10"^ 2 x 10"^ 6 X 10"'', 2 x 10"" 
from left to right. The asymptotes Dy yt'^, using the same 
values of Dy as in Fig. [5T] for the respective shear rates, are 
shown as thin dashed lines. 



Fig. 26. The quantities 3Sx'^{t,t^)/{2j^t^) and 
Sy^{t,tiu)/{2jt), testing their approach of /3y in Eqs. H76p 
and H67p for i -^ cxd. Symbols denote the stationary, lines 
the transient curves. Open symbols and solid lines show the 
shear-, filled symbols and dashed lines the gradient direction. 
The Peo numbers are 2 x 10"^ (squares), 2 x 10"^ (circles), 
2 X 10~^ (triangles) and 2 x 10"'' (diamonds). For the two 
largest shear rates, collapse of shear and gradient directions 
for long times is visible. The inset shows the curves for the 
gradient direction magnified where convergence of /3y(7) to 
the value of the master-curve for Peo -^ (indicated by a 
horizontal black line) is observed. 




10"" 10'^ 10"^ 10"' 10° 10^ 10^ 10^ 10" 



approach the scaling constant /3y{'y ^- 0) w 1.4, indicated 
by the horizontal black line. We emphasize again that this 
number miiquely describes all possible MSDs (in 2D) for 
long times. 

In theory, the MSD for the shear direction contains 
a term of order 0{t^) (compare Eq. ([77])). Fig. [57] shows 
this MSD in simulations, after subtraction of the t'^-term 
and division by t. A term of order ©(t^), which would 
manifest itself in a linear increase of the curves at long 
times, however, cannot be resolved. 

In Fig. [55] we finally show the long time diffusion co- 
efficients for the y-direction (as de&ied in Eq. (|67p ). as a 
function of shear rate. For large shear rates, the diffusiv- 
ities for the different densities are very close together, a 
behavior which is known also from the macroscopic shear 
viscosities [55]. As shear rate gets smaller, the diffusivities 



Fig. 27. Investigation of the next to leading term in the MSD 
for the shear direction in simulations at 77 = 81 (glass). Peo 
numbers are 6 x 10"\ 2 x 10"\ 6 x 10"^ 2 x 10"^ 6 x 10"^, 2 x 
10"^, 6 X 10"'', 2 X 10"'' from left to right. Shown are transient 
(solid) and stationary (dashed) MSDs after subtraction of the 
leading term oc t^ and subsequent division by t. 



for the liquid densities finally approach a constairt value 
given by the diffusivity of the unsheared suspension. These 
decrease with density |55]. On the glassy side, we observe 
the approach of the scaling regime, where the diffusivities 
are linear in shear rate. Simulation and theory agree with 
respect to all these findings. Quantitatively, there is a fac- 
tor of roughly 55 between theory and simulation, see the 
discussion of Fig. [5T] 
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Fig. 28. Long time diffusion coefficient as a function of Peo 
for different densities. Solid lines show the theoretical data for 
£ = — 10~^,0, 10~^ from top to bottom. Theoretical data has 
been shifted by a factor of 55 along the a::-axis. A dashed bar 
shows the slope of 1. Dashed horizontal lines show the Peo — >■ 
asymptote from the quiescent system. 



9 Summary 

We discussed some of the characteristic features of tagged 
particle dynamics for glasses under shear. The transient 
tagged particle density correlator shows strong imprints 
from the shear stress after switch on. Directly after the 
stress-overshoot, the correlation functions decay quickly 
and are superexponential. Nevertheless, after the stress 
has relaxed to its final value, they drastically slow down, 
leading to the appearance of a direction-dependent final 
shoulder. Despite the strong anisotropy of the applied flow 
field, the correlation functions show rather small effects of 
anisotropy. The mean square displacement of the tagged 
particle shows an effect known at low density as Taylor 
dispersion, which at the glass transition appears in modi- 
fied form to obey the scaling with shear rate. The coupling 
of the MSDs for shear and gradient directions is identical 
to the low density case. 

The extension of mode coupling theory to sheared sys- 
tems (MCT-ITT) allows to study the properties of the 
tagged particle correlator and the MSD. It captures many 
nontrivial effects (e.g. anisotropy of correlation functions, 
superexponential behavior, emergence of shoulders, scal- 
ing behavior in the glass for both stationary and tran- 
sient functions, Taylor dispersion), and gives quantitative 
predictions without adjustable parameters, where the re- 
sulting timescales are captured correctly within roughly 
one order of magnitude. Wc attribute the deviations in 
timescales to an underestimation of the stress-overshoot 
scenario in theory, as both correlators and MSDs do not 
speed up as strongly as the simulation curves after the 
stress maximum is passed. 
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A a-scaling equation 

To derive the equation for ^q+(t). wc start with the con- 
volution integral in Eq. PO)) . 



m^(i,i')a,'^^(i')- 



(89) 



In Eq. (|42p . we saw explicitly how the memory function 
depends on the two different times, namely by 



m^(t,i') =m^(^t,)(t-t'). 



(90) 



It depends only on accumulated strain ■^t' rather than on 
t' . It does so via the advected wavevector. That means 
that the dependence on t' is already a-scaling-like. Using 
this, we can rewrite the integral in Eq. (15^1) to 



VaieH0"''>(^*')(*-O)n(O. (91) 



la 9t' 

With this, the a scaling equation, Eq. ((50)) . follows. 

B Long time solution of the MSDs 

It is useful to rewrite the memory function m^{t, t') into a 
product of the part which depends only on the difference 
t ~ t' and the part which depends explicitly on t and t' , 

(92) 
where the function 

F(k,q(t'),t-t') = 



N 



-n^c' 



k(t-t') 



45fc<?^_q(to(i-i')<^k(i-i') (93) 



depends still explicitly on t' via the wavevector q(t'). How- 
ever, this dependence will vanish in the low q limit as used 
for the calculation of the MSDs. 



B.l Flow direction 

In order to find the long time solution of Eq. (|71l) we write 



lim 5x^{t) = at^ + hi^ + ct + 



(94) 



The form (|M1) can be justified knowing that the function 
F in Eq. ^^ decays to zero as i — i' -^ oo. A term f^ 
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(or higher powers in t, or fractional powers) does not ex- 
ist since the initial decay rate (|73[) does not contain such 
a term. The long time behavior is hence governed by the 
initial decay rate, a fact which is interesting because the 
long time behavior of the correlator ^q(i) is independent 
of the initial decay rate as 7 — )■ [29] . This is because the 
limits of i — > (X) and q — > do not commute. 
We first determine the coefficient a. For this, the leading 
long time (large t' and t") terms in the integral in Eq. (|71l) 
are needed. They are independent of the coefficient b. The 
equation for 6, on the other hand, will contain the coeffi- 
cient a and can hence only be solved afterwards. The lead- 
ing term of the first bracket in m^{t, t') is —'ytky{t—t'), the 
leading term of the fraction is given by —ky/{'^t'). With 
this, we get 



3at^ + I dt'^tm^it - t')£l^ = 272^2. (95) 



7^' 



We note that mPAt) appears. This equation can be treated 
with the following formula for Laplace transforms [43], 



B.2 Cross correlation 

In order to find the long time solution of Eq. (j79p . we have 
to find the long time behavior of 



1 - <?•", 

lim 2(2^ 



,w 



q^ 



Its equation of motion is given by. 



5xv{t)/2. 



(101) 



dt5xy + f dt'Y] [fc, -t- (1 - 'yt)kyit - t')] 
Jo k 



1 + (1 - it'y 
it) 



r, 



qe^+qe. 



q- 



(102) 



Again, we can only solve this equation for the long time 
contributions after making the following ansatz. 



lim Sxy{t) = a't^ + b't^ + c't + 



(103) 



C{tf{i)}{z)^--C{f{t)}{2 
oz 



(96) 



Using it, we find that the integral in Eq. (j95p contain also 
one term of order t^ (because dzm^iz) is finite as z — !• 0), 
which does not contribute to a. We find 
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3 4-3mO(z = 0)' 



(97) 



We must also consider the next order leading term as it 
will be needed in order to calculate the xy cross correla- 
tion. For the equation for the coefficient &, all long time 
terms proportional to t^ have to be collected, (note that 
some of the possible contributions vanish in the sum over 
k due to symmetries 13) , 



Wc note that the leading long time term in Eq. (J102I) is 
equal to the long time term of fe^, i.e., 

a = a. (104) 

Additionally to the terms in Eq. (|M|) . the equation for b' 
contains one contribution from the initial decay rate. The 
additional terms in Eq. (jl02p that come from the memory 
function exactly cancel each other in this order. We find 
for 6', 

^y T- (105) 

l + m"(z = 0) ^ ' 



b' = b- 



It is important that a' = a, leading to the cancellation of 
the i^-terms in Eq. ([71]). Putting the result of Eq. (fTU5|) 
into Eq. ^ leads to Eq. pi)) . 



m" (z = 0) 
2b + 26m°(z = 0) - 'Sad^mliz = 0) - 3a— ^^. = 0. 

(98) 
The promised dependence of 6 on a appears. Also, the 
off-diagonal memory function enters. 



ly{t)=Y.^y{t)k,F{\^,0,t). 



(99) 



We find for b, 



3a 



[9.m°(z = 0)]+<4f^ 



1 + toO(z = 0) 



(100) 



^ Although the memory function including the correlators is 
not isotropic under shear, it is still symmetric with respect to 
the origin, k = 0, since the system is symmetric with respect 
to the origin. 
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